3 dimensional matching / 3х-мерное сочетание

import random
from IPython.core.display import SVG, HTML
import pyomo.environ as pyo
from pysat.solvers import Solver
from pysat.formula import CNF 
import py_svg_combinatorics as psc
from ipywidgets import widgets, HBox
from collections import Counter
from pprint import pprint
from random import randint
import numpy as np
from IPython.display import IFrame
import IPython
from copy import copy
import os
import re
from pathlib import Path
nbname = ''
try:
    nbname = __vsc_ipynb_file__
except:
    if 'COCALC_JUPYTER_FILENAME' in os.environ:
        nbname = os.environ['COCALC_JUPYTER_FILENAME']    
title_ = Path(nbname).stem.replace('-', '_').title()    
IFrame(f'https://discopal.ispras.ru/index.php?title=Hardprob/{title_}&useskin=cleanmonobook', width=1280, height=300)

matchings3d = [
        [0, 'b', 77],
        [1, 'a', 88],
        [1, 'b', 12],
        [1, 'c', 12],
        [2, 'a', 12],
        [3, 'c', 77],
        [2, 'b', 12]
]

svg = psc.matchings2svg(matchings3d)
SVG(data=svg)

def print_solution(m):
    for v in m.component_data_objects(pyo.Var):
        if v.value and v.value > 0:
            print(str(v), v.value)

def get_model(matchings3d):
    m = pyo.ConcreteModel()
    m.E = [ (e[0], e[1], e[2]) for e in matchings3d ]

    m.V = sorted(set([(j, item) for sublist in matchings3d for j, item in enumerate(sublist)]))
    max_len = max([len(sublist) for sublist in matchings3d]) #3

    # Выбираем ли мы в оптимальное сочетание это гиперребро
    m.x = pyo.Var(m.E, domain=pyo.Binary)
    m.match_count = pyo.Objective(expr = sum( m.x[e] for e in m.E), sense=pyo.maximize)

    @m.Constraint(m.V)
    def нет_пересечениям(m, j, it):
        return sum(m.x[e] for e in m.E if e[j] == it) <= 1

    return m

m = get_model(matchings3d)
m.V, m.E

m.pprint()

solver = pyo.SolverFactory('cbc')
solver.solve(m).write()
print_solution(m)

print_solution(m)

selected = [j for j, e in enumerate(m.E) if m.x[e].value > 0]
selected

svg = psc.matchings2svg(matchings3d, selected)
SVG(data=svg)

cnf3 = CNF(from_clauses=psc.rand3cnf(3, 6))

from pysat.solvers import Solver

solver = Solver(bootstrap_with=cnf3)
res = solver.solve()
res 

print(solver.get_model())

cnf = CNF(from_clauses=psc.rand3cnf(2, 6))

cnf.clauses

def lit2vertex(cnf, lit, num=None):
    '''
       Отражение литерала в вершину 3д-сопоставления.
       Формула-представление для положительного или отрицательного литерала.
    '''
    n = cnf.nv
    x = abs(lit)
    v = f'x_{{{x}}}'
    if num is not None:
        v = f'x^{{{num}}}_{{{x}}}'
    if lit < 0:
        v = '¬' + v
    return v

def clause23dm(cnf, clause_num, with_garbage_collector=True):
    '''
       Скобка КНФ в набор 3d-сочетаний 
    '''
    match = []
    vn2count = Counter()  # счетчик встреч переменной в литералах
    for j, clause in enumerate(cnf.clauses): 
        clause_len = len(clause)
        for lit  in clause:
            vn = abs(lit)      # номер переменной
            vn2count[vn] += 1  # подсчитываем номер переменной по счету.
            if j == clause_num:
                b = lit2vertex(cnf, lit, vn2count[vn])
                match.append([b, f'y_{{c{j}}}', f'z_{{c{j}}}'])
                if with_garbage_collector:
                    for gn in range(clause_len-1):
                        nb = lit2vertex(cnf, -lit, vn2count[vn])
                        match.append([b,  f'y^{{g{j}}}_{{{gn}}}', f'z^{{g{j}}}_{{{gn}}}'])
                        match.append([nb, f'y^{{g{j}}}_{{{gn}}}', f'z^{{g{j}}}_{{{gn}}}'])
    return match

cnf.clauses

SVG(psc.rows2svgtable(clause23dm(cnf, 1, with_garbage_collector=False)))

svg = psc.matchings2svg(clause23dm(cnf, 0, with_garbage_collector=False))
SVG(data=svg)

match4clause = clause23dm(cnf3, 1)
SVG(svg1:= psc.matchings2svg(match4clause))

SVG(psc.rows2svgtable(match4clause))

# Не обращать внимание, отладочная ячейка с проблемами, пока не получается две SVG-картинки параллельно вывести.
#no_wrap_div = '<div style="white-space: nowrap">{}{}</div>'
#
#display(HTML(no_wrap_div.format(svg1, svg2)))

def variable23dm(cnf, var_num):
    '''
    Делаем гаджет колесо для заданной переменной.
    '''
    match = []
    k=0
    for clause in cnf.clauses:
        for lit  in clause:
            if abs(lit)==var_num:
                k += 1    
    n=k
    k=0
    for clause in cnf.clauses:
        for lit  in clause:
            if abs(lit)==var_num:
                vn = abs(lit)
                k += 1    
                bp = lit2vertix(cnf, vn, k)
                bn = lit2vertix(cnf, -vn, k)
                match.append([bp, f'y_{{{vn}}}^{{{k}}}', f'z_{{{vn}}}^{{{k}}}'])
                match.append([bn, f'y_{{{vn}}}^{{{(k % n)+1}}}', f'z_{{{vn}}}^{{{k}}}'])
    return match    

cnf.clauses 

m_ = variable23dm(cnf, 2)
SVG(psc.matchings2svg(m_))

def sat323dm(cnf):
    m = len(cnf.clauses)
    n = cnf.nv
    matchings = []
    for v in range(1, n+1):
        matchings += variable23dm(cnf, v)
    for j in range(m):
        matchings += clause23dm(cnf, j)
    return matchings

cnf.clauses

matchings = sat323dm(cnf)
svg = psc.rows2svgtable(matchings)
SVG(data=svg)

m = get_model(matchings)
solver.solve(m).write()
print_solution(m)

R = [e for j, e in enumerate(m.E) if m.x[e].value > 0]

R = [e for j, e in enumerate(m.E) if m.x[e].value > 0]
S = [j for j, e in enumerate(m.E) if m.x[e].value > 0]

#svg = psc.matchings2svg(matchings, selected=R)
SVG(psc.rows2svgtable(m.E, selected_rows=S))

S = [j for j, e in enumerate(m.E) if m.x[e].value > 0]
svg = psc.matchings2svg(matchings, selected=S)
SVG(data=svg)

svg = psc.rows2svgtable(matchings, selected_rows=S)
SVG(data=svg)

matchings

def matchsol2cnfsol(m):
    '''
    Конвертируем ЦЛП-решение для matchings в присваивания для CNF
    '''
    import re
    cnfsol = {}
    re_ = re.compile(r"(?P<neg>¬)?x\^\{\d+\}\_\{(?P<varn>\d+)\}")
    
    for e in m.E:
        if m.x[e].value > 0:
            if 'c' in e[2]:
                rm = re_.match(e[0])
                assert rm, e[0]
                varn = int(rm.group("varn"))
                val = 1
                if rm.group("neg"):
                    val = 0
                assert varn not in cnfsol or val == cnfsol[varn], (e, cnfsol[varn])
                cnfsol[varn] = val
    return cnfsol        

cnf.clauses

matchings = sat323dm(cnf)
m = get_model(matchings)
solver.solve(m).write()
print_solution(m)

SVG(psc.matchings2svg(variable23dm(cnf, 1)))

print_solution(m)

matchsol2cnfsol(m)

def solve_cnf_by_3dm(cnf):
    matchings = sat323dm(cnf)
    m = get_model(matchings)
    solver.solve(m)   
    sol = matchsol2cnfsol(m)
    return sol    

cnf = CNF(from_clauses=psc.rand3cnf(16, 9))
cnf.clauses, solve_cnf_by_3dm(cnf)